function Q = cubic_symm(q)
%Q = cubic_symm(q) -- Replicate q with the cubic symmetry operator.

S =   [[0    0    0    1];  ...
       [1    0    0    0];  ...
       [0    1    0    0];  ...
       [0    0    1    0];  ...
       [1    0    0    1] / sqrt(2);  ...
       [0    1    0    1] / sqrt(2);  ...
       [0    0    1    1] / sqrt(2);  ...
       [-1   0    0    1] / sqrt(2);  ...
       [0   -1    0    1] / sqrt(2);  ...
       [0    0   -1    1] / sqrt(2);  ...
       [1    1    0    0] / sqrt(2);  ...
       [-1   1    0    0] / sqrt(2);  ...
       [0    1    1    0] / sqrt(2);  ...
       [0   -1    1    0] / sqrt(2);  ...
       [1    0    1    0] / sqrt(2);  ...
       [-1   0    1    0] / sqrt(2);  ...
       [0.5  0.5  0.5  0.5];  ...
       [-0.5 -0.5 -0.5  0.5]; ...
       [0.5 -0.5  0.5  0.5];  ...
       [-0.5  0.5 -0.5  0.5]; ...
       [-0.5  0.5  0.5  0.5]; ...
       [0.5 -0.5 -0.5  0.5];  ...
       [-0.5 -0.5  0.5  0.5]; ...
       [0.5  0.5 -0.5  0.5]];

S = S(:, [4 1 2 3]);  % switch from mtex to libbingham format

   
Q = zeros(size(S));
for i=1:size(S,1)
    Q(i,:) = quaternion_mult(q, S(i,:));
end
